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fj , Abstract. An efficient Monte Carlo algorithm for the simulation of spin models 

^ ' with long-range interactions is discussed. Its central feature is that the number of 

operations required to flip a spin is independent of the number of interactions be- 
tween this spin and the other spins in the system. In addition, critical slowing down 
is strongly suppressed. In order to illustrate the range of applicability of the algo- 
rithm, two specific examples are presented. First, some aspects of the Kosterlitz- 
Thouless transition in the one-dimensional Ising chain with inverse-square inter- 
actions are calculated. Secondly, the crossover from Ising-like to classical critical 
e, behavior in two-dimensional systems is studied for several different interaction pro- 

files. 

O ' 1 Introduction 

For a long time, systems with long-range interactions have posed a great 
challenge to the application of Monte Carlo (MC) methods. Due to the large 
number of interactions that have to be taken into account, simulations are 
very time consuming and hence restricted to relatively small system sizes. 

^SJ , Alternatively, the interaction is truncated beyond a certain distance, which 

^O ' may introduce severe errors in the calculation. Thus, it is desirable to have 

an algorithm at one's disposal that allows the efhcient simulation of systems 
with long-range interactions without making any approximation except for 
the inherent statistical errors. A few years ago, an algorithm fulfilling these 
C^ ' requirements was indeed introduced for the simulation of spin models p| . It 

is based on the cluster method introduced by Swcndscn and Wang [p|| and 
has the crucial property that the number of operations required to add a 

' ^ I single spin to a cluster is independent of the number of interactions of this 

h^ ' spin with other spins in the system. Thus, the numerical effort to simulate 

(\ , ferromagnetic long-range models becomes comparable to that for short-range 

models. Close to the critical point the application of the algorithm becomes 
particularly profitable, since there one also benefits from the inherent reduc- 

j^ ' tion in critical slowing down in cluster algorithms. It is the purpose of the 

JH I present paper to discuss this algorithm in some more detail and to illustrate 

its power by means of some selected applications. 






o 



2 Erik Luijten 

2 Description of the Algorithm 

2.1 Conventional Wolff Cluster Algorithm 

To set the stage, we first briefly discuss the standard cluster algorithm for a 
system with short-range interactions. For simplicity, we will employ a formu- 
lation in terms of the Wolff cluster algorithm j^ rather than that of Swendsen 
and Wang (SW). It should be noted, however, that the modifications required 
for the simulation of long-range interactions pertain to the cluster-formation 
process only, and hence can equally well be combined with the SW variant. 
Since a detailed derivation of the Wolff algorithm can be found in numerous 
references, we restrict ourselves here to a schematic outline. Consider a d- 
dimensional lattice structure, where on each lattice site i a spin Si is placed, 
which can take the values ±1 and which has a ferromagnetic coupling K 
with its nearest neighbors. Thus, the system is described by the following 
Hamiltonian, 

f3H = -KJ2s^SJ, (1) 

where the sum runs over all pairs of nearest neighbors. At each Monte Carlo 
step, a cluster of spins with identical sign is flipped. The nonlocal character 
of this operation leads to a rapid decay of correlations. The crucial step is 
the identification of a cluster: 

1. Randomly choose a spin Sj from the lattice. 

2. Consider each spin that interacts with spin Si. It is added to the cluster 
with a probability p = I — exp(—2K), provided that it has the same sign 
as Si . 

3. Repeat step 2 in turn for each spin that is newly added to the cluster, 
where one now considers all spins that interact with this spin, rather than 
with Si . This is iterated until all neighbors of all spins in the cluster have 
been considered for inclusion. 



This cluster-building process is based on the Kasteleyn-Fortuin mapping | 
of the Potts model on a bond-percolation model. This mapping allows a 
generalization of the outlined procedure to systems with different interaction 
strengths Kj. In the second step, a spin is simply added with a probability 
Pj = 1 — exp{—2Kj) if the interaction strength with the spin Si is equal 
to Kj . As a side note, it is remarked that this implies that a cluster now not 
necessarily can be visualized as a contiguous collection of spins. 

Although one can thus easily devise a cluster algorithm for systems with 
long-range interactions, its efficiency rapidly decreases with an increasing 
number of interactions. Indeed, the typical value for the coupling K in such 
a system will be relatively small and hence each single spin that is considered 
for inclusion in the cluster will be added only with a small probability. So, a 
lot of operations are required to add a single spin to the cluster. 
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2.2 Efficient Construction of Clusters 

In order to illustrate our approach, we take the example of a one-dimensional 
chain with a spin-spin interaction that decays as a power-law, Kij = fr~- ~'' . 
The cluster-building process starts with a spin on a randomly chosen site i 
and all other spins in the system are added to the cluster with a probability 
p(si, Sj) — SsiSjPij, where pij = 1 — exp[— 2/|i — jl"'-^^'^']. For each spin that 
is actually added to the cluster, its address is also placed on the stack. When 
all spins interacting with the first one have been considered, we read a new 
spin from the stack and repeat the process, until the stack is empty. The spin 
from which we are currently adding spins will be called the current spin. In 
order to avoid considering each single spin for inclusion in the cluster, we 
split up the probability p{si, Sj) into two parts, namely the Kronecker delta 
testing whether the spins Si and Sj have the same sign and the "provisional" 
probability pij. This enables us to define the concept of the cumulative prob- 
ability C{j), from which we can read off which spin is the next one to be 
provisionally added to the cluster. 



C(j)^£F(n) 



(2) 



n=l 



with 



P{n) 



n(i 



(3) 



Pj = l — eyi'p{—2Kj) is an abbreviation forpoj (and Kj = Kqj), i.e., we define 
the origin at the position of the current spin. P{n) is the probability that 
in the first step n — 1 spins are skipped and the nth spin is provisionally 
added. Now the next spin j that is provisionally added is determined by a 
(pseudo)random number g d [0,1): j — 1 spins are skipped if C{j — 1) < 
g < C{j). If the jth spin has indeed the same sign as the current spin then 
Sj is added to the cluster. Subsequently we skip again a number of spins 
before another spin at a distance k > j ]s provisionally added. Due to the 
requirement k > j we must shift the function P, 



Pjik) = 



n (1-^™) 

771— J + 1 



Pk , 



(4) 



and Eq. (m is simply a special case of Eq. (m . The appropriate cumulative 
probability is now given by a generalization of Eq. (H) , 



n=j + l 



(5) 
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By using the specific form of the probabihty pij one finds that this reduces 
to 

Q(fc) = l-exp|-2 ^ i^„ , (6) 

\ n=j+l J 

i.e., the probabihty that the next spin that wiU be added hes at a distance in 
the range [j + 1, k] is given by an expression that has the same form as the 
original probabihty, in which the couphng constant is replaced by the sum of 
all the couplings with the spins in this range. 

We now consider two possibilities of calculating the distance k from a 
given Cj(k). A straightforward approach is the construction of a look-up 
table. This means that we explicitly carry out the sum in (|6|) for a large 
number of distances fc, up to a certain cutoff M , and store the results in a 
table. Then, after drawing a random number, we can derive the corresponding 
spin distance from this table. In principle we need for each value of j another 
look-up table containing the Cj{k). This is hardly feasible and fortunately not 
necessary, as follows from a comparison of Eqs. ^ and (|^). Namely (assume 
k>j), 



C{k) = Co(fc) = C{3) + 



n(i-p^) 



Cj{k) 



= C{j) + [l^C{j)]C,{k) (7) 

or Cj{k) = [C{k) - C(j)]/[1 - C{j)]. So we can calculate Q(fc) directly 
from C{k). In practice one realizes this by using the distance j of the previous 
spin that was provisionally added to rescale the (new) random number g to 
g' e [C(j), 1); g' = C{j) + [1 — C{j)\g. Since we only consider ferromagnetic 
interactions, liaXj^aa C{j) exists and is smaller than 1, cf. Eq. (^. 

This method is very fast, since we have to calculate all cumulative proba- 
bilities only once, but it has two major drawbacks. First, we can accommodate 
only a limited number of spin distances in our look-up table and must there- 
fore devise some approximation scheme to handle the tail of the long-range 
interaction, which is essential for the critical behavior in the case of slowly 
decaying interactions (small a) . This issue is addressed below. Secondly, this 
method is impractical in more than one dimension, as the number of distances 
for which the cumulative probability has to be calculated rapidly increases 
with the dimensionality of the system (for a fixed cutoff) . 

For interactions that can be explicitly summed there exists a powerful 
alternative. In those cases, Eq. (0) can be solved for fc, yielding an expression 
for the spin distance in terms of Cj{k), i.e., in terms of the random num- 
ber g. Especially for the study of critical phenomena this is often a feasible 
approach, since in many cases one can modify an interaction such it can be 
explicitly integrated by only adding irrelevant terms that do not affect the 
universal critical properties. As an example, we consider again the power-law 



Monte Carlo Simulation of Spin Models with Long-Range Interactions 



interaction Kij = f\i — jI^'^^*^'. The corresponding sum appearing in the 
right-hand side of (o) is (for j = 0) the truncated Riemann zeta function, 
which cannot be expressed in closed form. However, upon approximation of 
this sum by the integral 



E^« 



k 



f 



f 



k+4 



dx X 



-(l+a) 



(8) 



dx X ^^ 



(9) 



one still has an exact Monte Carlo scheme, but for an interaction 

/■K-jl + i 
K{\i^]\) = f 

Both interactions belong to the same universality class, so that, e.g., for the 
determination of critical exponents one is free to make the most convenient 
choice. Of course, all nonuniversal quantities, such as the critical temperature, 
will have different values. For the modified interaction, Eq. (0) reduces to 

where it should be noted that this is only the probability for adding a spin 
that lies in a single direction. We have not written j + 5 and fc + 5 instead 
of J and A:, because we prefer to use continuous coordinates, which only in 
the last stage are rounded to a lattice site. Equating Cj(fc) to the random 
number g we find 



Cj{k) = 1 — exp 



(10) 



k^ 



r^ + 



2/ 



In(i-g) 



(11) 



Rescaling of the random number is no longer required: the lowest value, 5 = 0, 
leads to a provisionally added spin at the same distance as the previous one, 
k = j.li g = Cj{oo) = l — cxp[—{2f/a)j^"'] the next provisionally added spin 
lies at infinity and thus g e [Cj(oo), 1) yields no spin at all. Once the coordi- 
nates of the next provisionally added spin have been calculated by rounding 
to the nearest integer coordinate, the periodic boundary conditions are ap- 
plied to map this coordinate onto a lattice site. For the following provisionally 
added spin, j is set equal to k (not to the rounded distance!) and a new k is 
determined. If no spin has been added yet, j is set to ^, the lowest possible 
distance. For higher dimensionalities, d — 1 additional random numbers are 
required to determine the direction in which the next spin is added. 

The above-mentioned limited size M of the look-up table can now be 
coped with as well: beyond the spin distance M the sum in (||) is approx- 
imated by an integral. Thus, if the random number g lies in the interval 
\C{M), C(oo)), the spin distance k is determined from a modified version 
of (|ll| ) , where the lower part of the integral is replaced by an explicit sum 

1\-" /I *' 1 ''"'^" 

k^ (M+-) +a\—lii{l-g) + Y,' 



V 



,1+cr 



n=l 



(12) 
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2.3 Discussion 



Finally, we briefly illustrate the efflciency of the algorithm by means of an 
example. If the probability p were equal for each spin pair, one out of p~^ spins 
would be added to the cluster, and it would take 0{p~^) ^ ©(L") operations 
per spin to update a configuration, compared to 0{L'^) operations per spin for 
a Metropolis algorithm. Taking into account the decrease in critical slowing 
down, we see that the efficiency of this method is typically a factor 0{L'^^^) 
larger than the conventional Monte Carlo algorithm. In the case of a large 
but finite interaction range R, such as for the two-dimensional equivalent- 



neighbor model discussed in Sec. 3.2, one gains a factor 0{R'^L'^), provided 



that critical slowing down is completely suppressed in the cluster algorithm. 
Far from the critical temperature, the factor L^ disappears, but one still has 
the advantage of a speed increase proportional to the number of interactions 
per spin. 

3 Applications 

3.1 The Inverse-Square Ising Chain 

As a first application, we consider the Ising chain introduced before, with 
spin-spin interactions that decay as fr~- ~'^ . Unlike the Ising chain with 
short-range interactions, in which no long-range order is possible for nonzero 
temperatures, this system exhibits a remarkably rich behavior. The inter- 
ested reader is referred to Ref. |6| for a review. In summary, the system 
exhibits mean-field-like critical behavior for < u < ^ and nonclassical 
critical behavior for ^ < a < 1, with critical exponents that depend on a. 
For (T > 1 the interactions are essentially short-ranged. At tr = 1, the so- 
called inverse-square Ising model, we have a very interesting situation: the 
spin chain displays a phase transition which is the one-dimensional analog of 
the Kosterlitz-Thouless (KT) transition in the two-dimensional XY model. 
It has close connections to a variety of physical applications, such as the 
Kondo problem, quantum tunneling in a two-state system coupled to a dis- 
sipative environment, quark confinement, etc. Although the KT transition 
has received an enormous amount of attention over the past decades, there 
are still open questions. Two peculiar properties of this transition are: (1) at 
the critical temperature T^ the order parameter exhibits a singular behavior 
superposed on a jump like one finds for a first-order transition; (2) the cor- 
relation length ^ and the susceptibility x diverge exponentially for T [ T^. 
e = Co exp[Bj/(T - T^Y] and x = Xo exp[B^/(r - T^f], with y ^ \. These 
and related critical properties turn out to be very difficult to verify in nu- 
merical simulations of two-dimensional models, because the finite-size effects 
decay only logarithmically. The onc-dimcnsional model, which now can be 
simulated with comparable efficiency, clearly offers a great advantage: rather 
than a linear system size of 0(10^), like for d = 2, one can reach system sizes 
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L = 0(10^) and thus approach the critical point much closer. Since a detailed 
description of the critical properties and their numerical determination is be- 
yond the scope of this paper, we restrict ourselves here to a study of the 
quantity ^ = Km?, where K is the coupling constant and m the magnetiza- 
tion density, and its application for the determination of Tc (or, equivalently, 
its inverse, the critical coupling Kc). It can be shown that ^ is the analog of 
the spin-wave stiffness in the 2D XY model. Just like the latter quantity is 
predicted to have a universal jump 2/ti at criticality, 'P is expected to have 
a jump of size ^ 0]. Indeed, in Figure |l|, where ^(if, L) is shown for system 
sizes up to L = 400 000, one can already clearly observe how such a jump 
develops with increasing system size. The superposed square-root singularity. 



^{K,oo) = W{K^, oo) -I- Cy/K -Kc + 0{K - Kc) 
is shown in the inset. 



(13) 
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Fig. 1. The quantity ^ = Km as a function of the couphng K, for system sizes 
10 < i < 400 000. The inset shows limz,^^ >P{K, L) for K > K^. 



We consider now three distinct ways to use this quantity for the deter- 
mination of Kc, which for phase transitions in this universality class is no- 
toriously difBcult. First, one may use the predicted singular behavior of ^ 
for K > Kc- To this end, the finite-size data at fixed couplings first must be 
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extrapolated to the thermodynaniic limit. By integrating the RG equations 
one finds M that ^{K, L) obeys a finite-size expansion of the form 



1'{K, L) = ^{K, oo) {l + aiL-2[^-i] + a^i""''^-'! + • • •} 



(14) 



where ^ = ^{K, oo)/^{Kc, oo) and the ellipsis denotes higher-order terms. 
The resulting estimates for ^{K,(X)) are then fitted to Eq. (^3|), in which 
<Z/(Xc,L) is kept fixed at i. This yielded K^ = 0.6552 (2). Secondly, <l'iK,L) 
in the high-temperature regime K < K^ may be used to estimate K,,. Indeed, 
within the finite-size regime one expects two types of corrections to scaling: 
corrections due to irrelevant fields, which decay as powers of 1/lnL, and 
temperature-dependent corrections which can be expanded in terms of L/S,. 
A least-squares fit of the numerical data to an expansion of the form 



>F{K,L) = <F{K,,oo) 



L 



02 I 7- 



bi 

InL 



(lnL)2 



(15) 



has yielded ^{K^.oo) = 0.496 (3) and K^ = 0.6548 (14). Fixing ^{K^,oo) 
at the predicted value ^ we found Kc = 0.6555 (4). Finally, a very straight- 
forward but remarkably effective approach is to fit a set of finite-size data 
for W{K,L) at fixed coupling to an expression of the form of Eq. (|l^), in 
which all temperature-dependent terms have been omitted, i.e., one assumes 
that K = Kc- Although in principle such a fit should only work at the true 
critical coupling, it turns out that least-squares fits of a good quality can 
be obtained over a range of couplings. However, the resulting estimate of 
\I'{K, oo) is a monotonously increasing function of K and Kc can be deter- 
mined from the requirement that 'l'{Kc,oo) = |, yielding K,. — 0.65515 (20). 



Table 1. Some estimates (in chronological order) for the critical coupling Kc of 
the inverse-squaxe Ising chain. 



Reference Kc 



This 



M 

PI 

[Lil 

[121 
[131 
[14| 
[151 
[161 
[171 
[181 

[[9 1 
work 



^ 0.562 

0.5 

0.590 (5) 

0.615 

> 0.61128 



Method 



~ 0.612 Exact summation -I- Pade approximants 

^ 0.635 RG 

~ 0.490 Extended mean-field approach 

0.657 (3) Series expansion 

> 0.441 Extended mean-field approach 

Analytical 

Coherent-anomaly method 

Variational method 

Cycle expansion 

Transfer matrix 

Extended mean-field approach 
6/tv'^ ~ 0.6079 RG (conjectured to be exact) 
0.6552 (2) MC 
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It is rewarding that the three different methods yield consistent estimates for 
the critical coupling. Table nl summarizes several estimates for Kc^ obtained 
by a variety of methods. Remarkably, most numerical results (of which only 
few carry an explicit estimate of the uncertainty) suggest a critical coupling 
around 0.61. Our estimate lies considerably higher and only agrees with the 
series-expansion result of Rcf. 111]. All methods that rely on an extension 
of mean-field theory clearly suffer from the extremely slow convergence to 
the thermodynamic limit. The exact conjecture of Rcf. [M (which is, in- 
terestingly, precisely twice as large as the mean-field result S/tt^) has been 
refuted. 

3.2 Crossover from Ising-Like to Classical Critical Behavior 

In order to illustrate that also for interactions with a large but strictly fi- 
nite range the presented algorithm offers great advantages, we consider the 
so-called "equivalent-neighbor" model. This is a simple generalization of the 
Ising model, in which each spin interacts equally strongly with all its neigh- 
bors within a distance Rm ■ In the limit Rm — > oo this model becomes equiva- 
lent to the exactly-solved mean-field model, whereas all systems with a finite 
interaction range belong to the Ising universality class. Since an exact solution 
for the latter case is lacking for d = 3, it is interesting to study the variation 
of critical properties as a function of i?,„ (cf. Ref. (20|). Another, experi- 
mentally very relevant, application of this model will be discussed here. As is 
well known, many thermodynamic properties show a characteristic power-law 
divergence upon approach of a critical point. These powers, or critical expo- 
nents, have universal values that only depend on a small number of global 
properties of the system under consideration. For example, binary mixtures, 
simple fiuids, and uniaxial ferromagnets all exhibit the same set of critical 
exponents as the three-dimensional Ising model. However, the corresponding 
power-law behavior is only observed asymptotically close to the critical point, 
whereas at temperatures farther away from T^ (but still relatively close to it) 
one may observe classical or mean- field- like critical behavior. This crossover 
can be explained in terms of competing fixed points of a renormalization- 
group transformation and is in principle well understood. However, unlike 
the critical exponents, for which accurate results have been obtained from 
series expansions, renormalization-group calculations, experiments, and nu- 
merical calculations, the precise nature of the crossover from one universality 
class to another is still a point of discussion. In particular, it is an unset- 
tled question to what extent this crossover is universal. There exist several 
field-theoretic calculations, but it has not yet been possible to verify their 
correctness by means of experiments. Measurements in the critical region are 
not only difficult, one also has to take great care to make the temperature 
distance to the critical point not too large, since one then would leave the 
critical region. As stated by the Ginzburg criterion, the crossover is a func- 
tion of i/G, where t = {T — Tc)/Tc is the reduced temperature and G a 
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system-dependent parameter. Throughout the crossover region, t has to be 
kept small, but t/G has to be varied over several decades. The large extent of 
the crossover region also emphasizes its experimental relevance: many mea- 
surements of critical exponents are actually made within this region rather 
than asymptotically close to the critical point and hence a detailed knowledge 
of crossover functions is required for a proper interpretation of the data. Since 
the Ginzburg parameter G is a function of the interaction range, it is well 
possible to construct such crossover functions from data obtained by means 
of simulations of systems with different interaction ranges, where for each 
system t is varied over a limited range only. In view of the large coordination 
numbers [0(10^)] that have to be reached within this approach, this is only 
feasible with an advanced algorithm. 

We concentrate here on one specific crossover function, namely that for 
the susceptibility in a two-dimensional (2D) system, both below and above 
the critical temperature (see Refs. Pllps] for a more detailed discussion of 
this topic). In the 2D Ising model, the susceptibility x diverges for i t as 
Af (— t)^^'"' and for t | as A^t^"^'^, where the amplitudes ^j are known 
exactly. Mean-field theory, on the other hand, predicts a susceptibility that 
for i t diverges as l/(— 2i) and for i | as 1/t. It is our aim to numerically 
determine the effective susceptibility exponent ^^^ = —d\nx/d\nt, which 
is expected to interpolate smoothly between 7/4 and 1. We have carried 
out MC simulations for square lattices with a maximum linear size L — 
1000 and interaction ranges up to i?~„ — 10000 (coordination number z = 
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Fig. 2. The effective susceptibility exponent 7^^ for T < Tc as a function of the 
crossover variable tR^ . 
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Fig. 3. The effective susceptibility exponent 7^ for T > Tc as a function of the 
crossover variable tR^. 



31416). It has been shown that the Ginzburg parameter G ex R^^, where 
the effective interaction range R^ — z~^J2j^i ki ~ i"jP (k* ~ i"jl ^ ^m) 
has been introduced to avoid lattice effects. For i < the susceptibility is 
calculated from the fluctuation relation x' = L'^{{m'^) ~ (|m|)^)/A;Br, whereas 
X = L'^{'W?)/k^T was used for i > 0. The exponent 7^ is obtained by 
numerical differentiation and shown as a function of t/G in Figs. |2| and ^, 
respectively. 

Several observations can be made concerning these graphs. In the first 
place, one notices that both functions smoothly interpolate between the Ising 
exponent and its classical counterpart. Also the expectation that the crossover 
takes several decades in the parameter t/G is confirmed. However, there is a 
striking difference between 7^^ and 7^^: whereas the latter exhibits a gradual 
decrease from 7/4 to 1 when the distance to the critical point is increased, the 
former drops much more rapidly and even goes through a minimum where 
7eff < 1, before reaching its asymptotic value. Interestingly, the occurrence 
of such a nonmonotonic variation has been inferred from RG calculations for 
three-dimensional systems. The results presented here are the first to show 
that, at least in the low-temperature regime, this behavior can actually be 
observed in systems as simple as the two-dimensional Ising model with an 
extended interaction range. 

The fact that the curves in Figs. and overlap for different values of R 
suggests that the crossover functions possess a certain degree of universal- 
ity. However, to what extent one may expect such universality in a region 
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where the correlation length no longer diverges is still an unanswered ques- 
tion. It has been suggested PSI that an additional length scale conies into 
play here. In order to investigate this possibility, we have carried out sim- 
ulations of two-dimensional models that are very similar to those studied 
above, except that the interaction profile has been modified p3. Instead of 
a constant interaction strength (block-shaped profile), a double-blocked po- 
tential was used, where each spin interacts with a strength Ki with its zi 
neighbors within a distance r < Ri (domain Di) and with a strength K2 
with its Z2 neighbors within a distance Ri < r < R2 (domain 02). In order 
to create a strong asymmetry, a strength ratio a = K1/K2 = 16 was cho- 
sen. The outer interaction range was kept fixed at R2 — -\/140, whereas the 
parameter Ri was varied. This allowed us to realize two additional, qualita- 
tively different situations: both Rx = \/A and i?i = -\/93 lead to an effective 
interaction range R « -x/SO [the previous definition for R is now generalized 
to i?^ = (aX^iGD ^"i + SieD ''f)/('^-^i + -^2)] , but in the former case the 
inner domain contains only 12 out of 436 interaction neighbors, compared to 
292 out of 436 in the latter case. This means that the integrated coupling 
ratio azi/ Z2, which indicates the relative contribution of the two domains to 
the total integrated coupling, is 0.45 in the first case and 32 in the second 
case. Extensive Monte Carlo simulations have been carried out for these and 
several other double-blocked interaction profiles. After a determination of the 
critical properties, we have calculated the crossover curve for the susceptibil- 
ity in the region t < for each of these profiles, in order to see whether the 
nonmonotonicity in Fig. is a peculiarity of the equivalent-neighbor model. 
The critical temperature, normalized by the critical temperature of the mean- 
field model, indeed displays a clear nonuniversality: for R\ = 93 the critical 
fluctuations turn out to be stronger suppressed than for R\ — 4, even though 
the effective interaction range is almost identical. The data for the suscepti- 
bility are shown in Fig. H. Since the critical amplitude of the susceptibility 
varies as R~^^^, a data collapse is obtained for x' /R^- One observes how all 
data perfectly coincide with the crossover curve for the equivalent-neighbor 
model. The deviations from the curve at the right-hand side of the graph 
are caused by the fact that, sufficiently close to Tc, the diverging correlation 
length is truncated by the finite system size. For the systems with Rf — 4 
and Ri — 93 (both with i?2 = 140) this happens in the figure at different 
temperatures, despite the fact that they have very similar values for the effec- 
tive range R. The reason for this is that the data points pertain to different 
system sizes, viz. L — 1000 and L — 300, respectively. The inset shows that 
for the same system size (L — 300) the data points for both systems virtu- 
ally coincide, even in the finite-size regime. The crossover function for the 
connected susceptibility appears to be insensitive to the introduction of an 
additional length scale in the interaction profile. We hence view this graph 
as a strong indication that crossover functions possess a considerable degree 
of universality. 
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Fig. 4. The crossover function for the connected susceptibility below Tc, for the two- 
dimensional Ising model with a double-blocked interaction profile, x = 'X /C{R), 
where C{R) is a range-dependent correction factor that accounts for the fact that 
the critical amplitude for small R deviates from the asymptotic range dependence. 
This only introduces a shift along the vertical axis. The solid curve indicates the 
crossover function for the block-shaped interaction profile and the dashed lines mark 
the mean-field ("MF") and the Ising asymptote. The numbers in the key refer to 
the values for R^ and R^ for each interaction profile. 



4 Outlook and Conclusion 



We have discussed a cluster algorithm for spin models with long-range inter- 
actions that is several orders of magnitude more efficient than conventional 
algorithms. Its usefulness has been demonstrated for a system with power-law 
interactions as well as for systems with medium-range interactions. In both 
cases, physical results have been obtained that until now could not be gener- 
ated within a reasonable amount of computing time. However, the scope of 
the algorithm goes far beyond what could be presented here. As far as purely 
ferromagnetic interactions are concerned, it can be efficiently applied to any 
interaction profile, including anisotropic ones. Also the order parameter must 
not necessarily be Ising-like: the extension to general 0{n) models m is simi- 
lar to that for the original Wolff algorithm and the first application to a Potts 



model has already been published |25 



The generalization to a situation in which additional antifcrromagnetic 
interactions are present is, in principle, also possible. Whereas competing 
interactions will move the system away from the percolation threshold (and 



14 Erik Luijten 

thus critical slowing down will no longer be optimally suppressed), one still 
has the advantage that not every individual spin has to be considered for 
inclusion in the cluster. Finally, it would be interesting to see whether this 
algorithm can be useful in other fields of physics where cluster algorithms 
have come to flourish. 
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